High order three part spht symplectic integrators: 
Apphcation to the disordered discrete nonhnear Schrodinger equation 
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While symplectic integration methods based on operator splitting are well established in many 
branches of science, high order methods for Hamiltonian systems that split in more than two parts 
have not yet been studied in detail. We demonstrate ways to construct high order symplectic 
integrators for Hamiltonian systems that can be split in exactly three integrable parts. Using these 
techniques for the integration of the disordered, discrete nonlinear Schrodinger equation, we show 
that three part split symplectic integrators are more efficient than other numerical methods for the 
long time integration of multidimensional systems, with respect to both accuracy and computational 
time. 

PACS numbers: 05.45.-a, 45.10.-b, 02.60. Jh, 05.60.Cd 



Following the time evolution of a dynamical system 
is generally accomplished by solving its corresponding 
equations of motion. If, for instance, the system under 
consideration can be described by an autonomous Hamil- 
tonian function H{q,p), virith q, p respectively being vec- 
tors of the generalized coordinates and momenta, the 
equations of motion can be readily derived via Hamilton's 
equations. One then attempts to determine the solution 
x{t) — {q{t),p{t)), t > 0, for any given initial condition 
x(0). Formally this solution can be described by the ac- 
tion of the operator e^^" , with Lh = X] i -^pi ^g, ~ ^q. dp- , 
on the initial condition, i.e. x{t) = e*"^" x{Q). The Hamil- 
tonian is said to be integrable if the action of this opera- 
tor is known explicitly and the solution of the Hamilton 
equations of motion can be written in a closed, analytic 
form. Unfortunately, this task is rarely possible, but in 
most cases the true solution can be approximated numer- 
ically. General purpose numerical integration methods 
for ordinary differential equations are capable of provid- 
ing such approximations. 

In this respect, the so-called symplectic integration 
techniques are of particular interest, as they are ex- 
plicitly designed for the integration of Hamiltonian sys- 
tems (see, for example. Chap. VI of [1| and refer- 
ences therein). Assume that H{q,p) can be written as 
H{q,p) = A{q,p) + B{q,p), so that the action of oper- 
ators e*^-^ and e*^^ is known, and the solution of their 
Hamilton equations of motion can be written analyti- 
cally, while e'^^" does not permit a closed analytical 
solution of its equations of motion. Then, a symplec- 
tic scheme for integrating the equations of motion from 
time t to time t-\-T consists of approximating the operator 
e'^^" = Q'^iLA+Ls) i^y product of j operators e"^'"^^^ and 
^diTLs ^ which represent exact integrations of Hamiltoni- 



ans A{q,p) and B{q,p) over times c^t and d^r respec- 
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constants Ci and di are appropriately chosen to increase 
the order of the remainder of this approximation. In 
practice, using this symplectic integrator (SI) we approx- 
imate the dynamics of the real Hamiltonian H — A + B 
by a new one, K ^ A + B + ©(t"), introducing an error 
term of order r" in each integration step, and the SI is 
said to be of order n. 

By their construction Sis preserve the symplectic na- 
ture of the Hamiltonian system and keep bounded the 
error of the computed value of H (which is an inte- 
gral of the system, commonly referred as the 'energy') 
irrespectively of the total integration time. Generally, 
this is not the case with non-symplectic integration al- 
gorithms. As a consequence. Sis permit the use of rela- 
tively large integration time steps r for acceptable levels 
of energy accuracy, resulting in lower CPU time require- 
ments. Recently, it was shown that Sis are also highly 
efficient in the integration of the variational equations 
needed for the computation of chaos indicators like the 
inaximum Lyapunov Characteristic Exponent (mLCE) 
;3i and the Smaller (SALI) and Generalized Alignment 
Index (GALI) Q when using the so-called 'Tangent Map' 
method Due to these benefits. Sis became a stan- 
dard technique in Hamiltonian dynamics with particular 
importance in long time integrations of multidimensional 
systems. Several Sis of different orders based on operator 
splitting have been developed over the years by various 
researchers 

In many cases of practical interest the Hamiltonian 
can be written as a sum of the system's kinetic energy 
T{p), dependent only on the momenta p, and the poten- 
tial V{q), dependent only on the positions q. In such 



2 



cases, the obvious choice for the apphcation of a SI is to 
consider A = T(p) and B = V(q). Yet in many phys- 
ical problems the corresponding Hamiltonian cannot be 
split in two integrable parts. The question then naturally 
arises whether it is possible to exploit the advantages of 
Sis for such systems as well. 

In this paper we present a systematic way to construct 
efficient, high order Sis for Hamiltonians that can be 
split in three integrable parts. Particular cases of sec- 
ond order three part split Sis, connected with astronom- 
ical problems, have been reported in literature llj. In 



these works, the considered Hamiltonians were expressed 
as if = A{q,p) + B{q,p) + C{q,p), the action of opera- 
tors e'^^'^, e"^^^ and e'^^'^ was analytically obtained, the 
second order SI of 5 steps 



'solution A' in [6[ 

ABC^{t) = ABC^{w3t)ABC^{w2t)ABC^{wit)x 
X ABC"^ {wot)ABC^ {wit)ABC^ {w2t)ABC^ (wgr) 

(4) 

having 29 steps. The exact values of Wi, i = 0, 1, 2, 3 can 
be found in [1, Chap. V, Eq. (3.11)] and @. 

In order to investigate the efficiency of the different SI 
schemes we choose a multidimensional Hamiltonian sys- 
tem describing a one-dimensional chain of coupled, non- 
linear oscillators. In particular we consider the Hamil- 
tonian of the disordered discrete nonlinear Schrodinger 
equation (DNLS) 

-Hd = E^'I^'I' + f l^'l' - (V^;+i^; +^r+i^O, (5) 
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was constructed, and its performance was studied. This 
integrator represents the simplest form of a symmetric SI 
that can be constructed for a Hamiltonian which splits in 
three distinct parts. Surprisingly enough, no systematic 
attempts for obtaining general high order three part split 
Sis have been performed to date, although integrator dU) 
showed a good numerical performance. Nevertheless, the 
idea of three part split was implemented for constructing 
second and fourth order integration schemes for a par- 



ticular complicated molecular model jl2| . while recently 
an adaptation of high order three part split methods was 
attempted to treat a specific astronomical problem Q . 

Our study fills this void by presenting a systematic 
way to construct higher order three part split Sis, based 
on the composition technique proposed by Yoshida 0. 
Starting from a SI ^^"(r) of order 2n, we can construct 
a SI S'2"+2(t) of order 2n + 2, as 



S 



2n+2 



(r) = S^'\z,t)S'^{zot)S'^{z,t), 



(2) 



with zo = -2i/(2"+i) /[2 - 2i/(2"+i)] and 

zi = l/[2 — 2^/('^"+^']. Applying this procedure to 
the second order SI ([Ij we obtain the fourth order SI of 
13 steps 

ABC^{t) ^ ABC^{xiT)ABC^{xaT)ABC^{xiT), (3) 

with xo = -2^/3/(2 - 2i/3) and xi = 1/(2 - 2^/^). 

Equation ^ can be used repeatedly to construct 
higher order three part split Sis. Although such a pro- 
cedure for constructing arbitrary Sis of even order with 
exact coefficients is straightforward, it is not optimal with 
respect to the number of required steps. As was already 
pointed out in Q, alternative methods can be applied 
to obtain more economical integrators of high order, al- 
though the new coefficients can no longer be given in 
analytical form. Several sixth order Sis of this kind were 
presented in [6] . Here, we consider one corresponding to 



with complex variables ipi , lattice site indices I and non- 
linearity strength /3 > 0. The random on-site energies e/ 
are chosen uniformly from the interval [— ^] , with W 
denoting the disorder strength. This model has two inte- 
grals of motion, as it conserves both the energy ([5]) and 
the norm S — J^i IV'/ ^^nd has been extensively investi- 
gated in order to determine the characteristics of energy 
spreading in disordered systems fl3l - [l7| . These studies 
showed that the second moment, 7712, of the norm dis- 
tribution grows subdiffusively in time t, as t", and the 
asymptotic value a = 1/3 of the exponent was theoreti- 
cally predicted and numerically verified. Currently open 
questions on the dynamics of disordered systems concern 
the possible halt of wave packet's spreading for i — > 00 
18 1 , as well as the characteristics of its chaotic behavior 
19|. Thus, providing the means to perform accurate long 
time simulations for the DNLS model within reasonable 
amounts of computational time is essential. 

Applying the canonical transformation 
i'l = in + ipi)/V2, "^i = {m - one can split 

([S]) into a sum of there integrable parts A, B and C as 
follows 



(6) 

where qi and pi are respectively generalized coordinates 
and momenta. For these three parts the propagation of 
initial conditions {qi,pi) at time i, to their final values 
{q'i,p[) at time t + t is given by the operators 



q'l = qi cos(a/r) + pi sin(a;T) 
p'l = pi cos(a;T) - qi sin(Q;;T) ' 



„tLb 



[pi-i +Pi+i)t 



Pi = PI 

Qi = n 

q'l = qi 

p'l = pi + {qi-i +qi+i)T 



(7) 



(8) 



(9) 
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with ai = ei + l3{qf +pf)/2. Thus, the DNLS model 
represents an ideal test case for our aforementioned three 
part split Sis. 

In order to evaluate the efficiency of the Sis (H]), 
and Q we compare their performance to that of other 
numerical techniques. 



In |14| - |l7l | numerical integration 
schemes based on traditional two part split Sis were ap- 
plied for the integration of Hamiltonian ([6]). These ap- 
proaches were based on the split of ([6]) in two parts as 
A = A and B — B + C, and the application of second or- 
der Sis of the so-called S AB A-family [l^ . In our study 
we implement the second order SI SABA2 using the split 
Hd — A + B. The integration of the A part is performed 
according to ([T]), while different approaches for approx- 
imating the action of e'^^'^ = e^^^+'^ are followed. In 



14i |15| a numerical scheme based on Fourier transforms 
was implemented (see appendix of 15| for more details) 
leading to a second order integrator with 5 steps, which 
we name SIFT^ in the following. Another approach is to 
split the B part in two integrable parts as S = i? -f C 
and use the SABA2 SI to approximate its solution. This 
means that we perform two successive two part splits in 
order to integrate Ho- This approach leads to a second 
order SI with 13 steps which we name SS^ (this scheme 
corresponds to the PQ method used in [l3])- In addition, 
based on two part split Sis we construct an integrator of 
order higher than two by applying the composition pro- 
cedure ([2]) to the SS^ integrator [i^l- In this way we 
construct a fourth order SI with 37 simple steps which, 
to the best of our knowledge, has never being used be- 
fore for the integration of the DNLS system, and we call 
it SS'^. 

Of course one can also use any general purpose non- 
symplectic integrator for the integration of ([6|) . One dis- 
advantage of such techniques is that different epochs of 
the system's evolution are computed with different accu- 
racy since these integrators do not keep the energy error 
bounded, but increase it as time increases. In particular 
for the DNLS model considered here the later stages of 
its evolution, which are of most importance since we are 
mainly interested in the asymptotic behavior of the sys- 
tem, are computed less accurately. As a representative of 
non~symplectic integrators we consider here the variable 
step Runge-Kutta method called DOP853, whose perfor- 
mance is controlled by the so-called one-step accuracy 5 

da. 

In order to compare the performance of the various 
integration schemes we consider a particular disorder re- 
alization of the DNLS model dH) with N = 1024 lattice 
sites. We fix the total norm of the system to S* = 1, 
and following [l^ we initially excite homogeneously 21 
central sites by attributing to each one of them the same 
constant norm, but with a random phase, while for all 
other sites we set qi{Q) — pi{Q) — 0. Due to the nonlin- 
ear nature of the model the norm distribution spreads, 
keeping of course the total norm S = ^i{qf + pf)/2 



constant (5* = 1). The performance of the integra- 
tion schemes is evaluated by their ability to (a) repro- 
duce correctly the dynamics, which is reflected in the 
subdiffusive increase of m2{t), (b) keep the values of 
the two integrals Hjy, S constant, as monitored by the 
evolution of the absolute relative errors of the energy 
Er{t) = liHoit) - i?D(0))/i/z?(0)|, and norm Srit) = 
\{S{t) - 5(0))/5'(0)|, and (c) reduce the required CPU 
time Tc{t) for the performed computations. 

Results obtained by the second order Sis ABC^, SS^ 
and SIFT^ and the non-symplectic integrator DOP853 
are presented in Fig. [1] These integration methods cor- 
rectly describe the system's dynamical evolution since for 
all of them the wave packet's ?7i2 shows practically the 
same behavior (Fig.[T^). The time steps r of the three Sis 
were chosen so that all of them keep the relative energy 
error practically constant at Er « 10~^ (Fig- [lb)- Since 
we are interested in the accurate long time integration 
of the DNLS model we use 5 = 10^^^ for the implemen- 
tation of the DOP853 integrator. For t « 10^ (which 
can be considered as a typical final integration time for 
long time simulations), this choice results practically in 
the same energy error obtained by all other tested in- 
tegrators. From Fig. [TJ; we see that the relative norm 
increases for all used methods, exhibiting larger 
values yet lower increase rates, for the ABC^ and SS^ 
Sis. Nevertheless, our results indicate that all methods 
can keep 5*^ to acceptable levels (e.g. Sr ^ 10~^), even 
for long time integrations. From Fig. [TJl we see that the 
SIFT^ integration scheme is the most efficient one with 
respect to the CPU time needed for obtaining the results 

of Fig. m 

For this reason we use the SIFT^ SI as a reference 
method, and compare in Fig. [2] its results with the ones 
obtained by our higher order Sis: the SS"*, ABC^ and 
ABC^ methods. These higher order Sis reproduce cor- 
rectly the evolution of m2 (Fig. [2^), keep Er « 10~^ 
(Fig. Wp), and show a slow increase of Sr with values 
remaining acceptably small (Fig. The ABC' and 

ABC^ Sis require less CPU times respective to the SIFT^ 
method (Fig. [5Jl), with ABC^ showing the best perfor- 
mance. From the results of Fig. [2] we see that using the 
ABC^ with T = 0.15 we need ~ 1.5 times less CPU time 
than the SIFT^ with r = 0.05. Although one might 
argue that this CPU time gain factor is not too big, 
we should keep in mind that long time simulations up 
to t = 10^ - 10^ of the DNLS model with N - 1000 
sites could require (depending on the particular com- 
puter used) up to ~ 10 days of computations. Thus a 
gain factor close to 2 is practically significant as it can 
considerably reduce the computation time. Our results 
indicate that the construction of efficient triple split Sis 
can allow the integration of the DNLS for longer times, 
and numerically tackle questions about the asymptotic 
behavior of wave packets. 

In summary, we presented ways to use Sis for Hamil- 
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1. Results for the integration of Hd ((6| by the second 
order Sis ABC^ for r = 0.005, SS^ for r = 0.02, SIFT^ for r = 
0.05 [(g) green; (b) blue; (bl) black], and the non-symplectic 
integrator DOP853 for 5 — 10^^^ [(r) red]: time evolution 
of the logarithm of (a) the second moment m2(t), (b) the 
absolute relative energy error Er{t), (c) the absolute relative 
norm error Sr{t), and (d) the required CPU time Tc{t) in 
seconds. 



tonian systems that do not split in two integrable parts, 
as traditional symplectic methods require, but in three. 
For such systems we constructed high order three part 
split Sis applying a systematic way for their creation: the 
composition method of @, and emphasized their practi- 
cal importance. In particular, we showed that such three 
part split Sis are more efficient numerical schemes than 
other symplectic and non-symplectic methods in terms 
of both accuracy and CPU time requirements, especially 
for the long time integration of multidimensional systems 
like the DNLS model. 

We hope that our results will draw the interest of the 
community in the construction of three part split Sis, 
and will initiate future research both for the theoretical 
development of new integrators of this type, as well as for 
their applications to different dynamical systems. Keep- 
ing in mind that such Sis can provide efficient numeri- 
cal schemes for the long time integration of Hamiltonian 
systems with many degrees of freedom (like the DNLS 
model), it would be interesting to investigate if the pos- 
sible addition of a corrector term can improve their accu- 
racy's done for traditional two part split methods (see 
e.g. 
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Max Planck Institute for the Physics of Complex Sys- 
tems in Dresden for its hospitality during his visits in 
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FIG. 2. Results for the integration of Hd (O by the second 
order SI SIFT^ for r = 0.05 [(bl) black], the fourth order Sis 
SS" for r = 0.1, ABC" for r = 0.05 [(b) blue; (g) green], and 
the sixth order SI ABC** for r = 0.15 [(r) red]. The panels 
are as in Fig. [1] Note that in panel (d) the black and blue 
curves practically overlap. 
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